Unify anisotropic viscosity Stokes assemblers / Fix anisotropic viscosity cookbook #6728
Conversation
gassmoeller
force-pushed
the
unify_anisotropic_viscosity_code
branch
from
7d29297 to
dbf80e1
Compare
|
/rebuild |
gassmoeller
force-pushed
the
unify_anisotropic_viscosity_code
branch
from
dbf80e1 to
e544695
Compare
|
/rebuild |
| subsection Function | ||
| set Variable names = x,z | ||
| set Function constants = pi=3.1415926; | ||
| set Function constants = pi=3.1415926 |
There was a problem hiding this comment.
Did that yield an error? If not, our parsing is not robust enough...
There was a problem hiding this comment.
Yes, the semicolon yields the following error (with deal.II v9.7.0-rc1). I guess that it worked with an older deal.II version when the cookbook was created (but have not tested specific versions yet). It is also interesting that the first error complains about the function expression, but the actual problem is inside the function constants.
ERROR: FunctionParser failed to parse
'Initial composition model.Function'
with expression
'0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02)); (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? sin(45*pi/180) : 0.0; (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? cos(45*pi/180) : 0.0;'
More information about the cause of the parse error
is shown below.
---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
to avoid a possible deadlock.
---------------------------------------------------------
----------------------------------------------------
Exception 'ExcMessage("Can't convert <" + s + "> to a double.")' on rank 0 on processing:
--------------------------------------------------------
An error occurred in line <678> of file </home/rene/software/deal.II-9.5.2/tmp/unpack/deal.II-v9.7.0-rc1/source/base/utilities.cc> in function
double dealii::Utilities::string_to_double(const std::string&)
Additional information:
Can't convert <3.1415926;> to a double.
Stacktrace:
-----------
#0 /home/rene/software/deal.II-9.5.2/deal.II-v9.7.0-rc1/lib/libdeal_II.g.so.9.7.0-rc1: dealii::Utilities::string_to_double(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&)
#1 /home/rene/software/deal.II-9.5.2/deal.II-v9.7.0-rc1/lib/libdeal_II.g.so.9.7.0-rc1: dealii::Functions::ParsedFunction<2>::parse_parameters(dealii::ParameterHandler&)
#2 ./aspect: aspect::InitialComposition::Function<2>::parse_parameters(dealii::ParameterHandler&)
#3 ./aspect: aspect::InitialComposition::Manager<2>::parse_parameters(dealii::ParameterHandler&)
#4 ./aspect: aspect::Simulator<2>::Simulator(ompi_communicator_t*, dealii::ParameterHandler&)
#5 ./aspect:
#6 ./aspect: main
--------------------------------------------------------
Aborting!
----------------------------------------------------
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
There was a problem hiding this comment.
A long time ago, we parsed numbers by calling a function that just tried to read a number from a string. That succeeds with 3.141;, it just doesn't eat the whole string. We only fixed that a couple of years ago.
There was a problem hiding this comment.
See dealii/dealii#19003. I now realize the problem is that the semicolon is not the correct separator for function constants (it should be a comma). I still think this probably used to work at some point, but does not anymore. Not sure if we want to do something about it, or leave it at the current behavior (since the documentation clearly states constants are separated by comma).
bangerth
left a comment
bangerth
left a comment
There was a problem hiding this comment.
I still have to look at the .cc file for the assembly, but here are a few comments already.
| * anisotropic viscosity tensor. If the material model provides a stress- | ||
| * strain director tensor, then the strain-rate is multiplied with this |
There was a problem hiding this comment.
| * anisotropic viscosity tensor. If the material model provides a stress- | |
| * strain director tensor, then the strain-rate is multiplied with this | |
| * anisotropic viscosity tensor. If the material model provides a stress-strain | |
| * director tensor, then the strain-rate is multiplied with this |
| const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity = | ||
| material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>(); | ||
|
|
||
| for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q) | ||
| { | ||
| // If there is an anisotropic viscosity, use it to compute the correct stress | ||
| const SymmetricTensor<2,dim> &directed_strain_rate = ((anisotropic_viscosity != nullptr) | ||
| ? | ||
| anisotropic_viscosity->stress_strain_directors[q] | ||
| * material_model_inputs.strain_rate[q] | ||
| : | ||
| material_model_inputs.strain_rate[q]); |
There was a problem hiding this comment.
Would it make sense to just always expect that a director is given when using this heating model?
| namespace aspect | ||
| { | ||
| namespace MaterialModel | ||
| { | ||
| /** | ||
| * Additional output fields for anisotropic viscosities to be added to | ||
| * the MaterialModel::MaterialModelOutputs structure and filled in the | ||
| * MaterialModel::Interface::evaluate() function. | ||
| */ | ||
| template <int dim> | ||
| class AnisotropicViscosity : public NamedAdditionalMaterialOutputs<dim> | ||
| { | ||
| public: | ||
| AnisotropicViscosity(const unsigned int n_points); |
There was a problem hiding this comment.
The assemblers directory is a bit of a funny place to find the declaration of a material model (or named additional output, as it may be). Have you thought moving it into the respective directory?
| namespace | ||
| { | ||
| template <int dim> | ||
| std::vector<std::string> make_AnisotropicViscosity_additional_outputs_names() | ||
| { | ||
| std::vector<std::string> names; | ||
|
|
||
| for (unsigned int i = 0; i < Tensor<4,dim>::n_independent_components ; ++i) | ||
| { | ||
| TableIndices<4> indices(Tensor<4,dim>::unrolled_to_component_indices(i)); | ||
| names.push_back("anisotropic_viscosity"+std::to_string(indices[0])+std::to_string(indices[1])+std::to_string(indices[2])+std::to_string(indices[3])); | ||
| } | ||
| return names; | ||
| } | ||
| } |
There was a problem hiding this comment.
Modules don't allow for anonymous namespaces in header files, and it's generally a bad idea anyway. I think that you're using this function only immediately below in one place, so it should be easy to move this code into its only caller.
| template <int dim> | ||
| AnisotropicViscosity<dim>::AnisotropicViscosity (const unsigned int n_points) | ||
| : | ||
| NamedAdditionalMaterialOutputs<dim>(make_AnisotropicViscosity_additional_outputs_names<dim>()), |
There was a problem hiding this comment.
Just convert the function above into a lambda function that you can declare and call in-place here.
| template <int dim> | ||
| std::vector<double> | ||
| AnisotropicViscosity<dim>::get_nth_output(const unsigned int idx) const | ||
| { | ||
| std::vector<double> output(stress_strain_directors.size()); | ||
| for (unsigned int i = 0; i < stress_strain_directors.size() ; ++i) | ||
| { | ||
| output[i]= stress_strain_directors[i][Tensor<4,dim>::unrolled_to_component_indices(idx)]; | ||
| } | ||
| return output; | ||
| } |
There was a problem hiding this comment.
This function could probably just as well live in a .cc file.
| namespace Assemblers | ||
| { | ||
| /** | ||
| * A class containing the functions to assemble the Stokes preconditioner. |
There was a problem hiding this comment.
| * A class containing the functions to assemble the Stokes preconditioner. | |
| * A class containing the functions to assemble the Stokes preconditioner for the case of anisotropic viscosities. |
| * A class containing the functions to assemble the compressible adjustment | ||
| * to the Stokes preconditioner. |
There was a problem hiding this comment.
| * A class containing the functions to assemble the compressible adjustment | |
| * to the Stokes preconditioner. | |
| * A class containing the functions to assemble the compressible adjustment | |
| * to the Stokes preconditioner for the case of anisotropic viscosities. |
| * This class assembles the terms for the matrix and right-hand-side of the incompressible | ||
| * Stokes equation for the current cell. |
There was a problem hiding this comment.
same here, and perhaps also below
| }; | ||
|
|
||
| /** | ||
| * This class assembles the term that arises in the viscosity term of Stokes matrix for |
There was a problem hiding this comment.
| * This class assembles the term that arises in the viscosity term of Stokes matrix for | |
| * This class assembles the term that arises in the viscosity term of the Stokes matrix for |
This is a copy-paste mistake -- the grammatically wrong code is already in line include/aspect/simulator/assemblers/stokes.h:81.
| std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity = | ||
| scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>(); |
There was a problem hiding this comment.
Could this variable be const?
| const double eta = scratch.material_model_outputs.viscosities[q]; | ||
| const double one_over_eta = 1. / eta; | ||
|
|
||
| const bool use_tensor = (anisotropic_viscosity != nullptr); |
There was a problem hiding this comment.
This variable should move out of the loop. But more generally, shouldn't this assemble only be called for the anisotropic case? We could assert that at the top, and then simplify a lot of the code below. (I recognize the history of the code as what we used to have, where we catered to both the anisotropic and isotropic case, but perhaps we don't need to carry this provenance around with us for forever.)
| template <int dim> | ||
| void | ||
| StokesCompressiblePreconditionerAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
There was a problem hiding this comment.
same comments for this function as for the one above.
| template <int dim> | ||
| void | ||
| StokesIncompressibleTermsAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
| template <int dim> | ||
| void | ||
| StokesCompressibleStrainRateViscosityTermAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
|
Thanks for the comments I will fix them. Originally I just left the code as it was to simplify review, but you are of course right that I may as well fix them immediately. |
gassmoeller
force-pushed
the
unify_anisotropic_viscosity_code
branch
from
a266d9c to
51c647a
Compare
|
I think I addressed all your comments except for:
I moved the relevant code out of the header into a source file (that is where it should be anyways). I would know how to convert this into a lambda function, but I feel that would reduce the readability of the code a lot. The return value of this function is used as an argument in a call to the constructor of a base class. Introducing and immediately calling a lambda function here would look very complicated, and at least I would prefer the current solution. Or did I misunderstand your suggestion? If the anonymous namespace is the problem, I would be happy to rename the namespace, e.g. to Everything else is addressed, in particular I moved the Material Model Additional Outputs object into a new directory |
bangerth
left a comment
bangerth
left a comment
There was a problem hiding this comment.
This looks good. The issue with the lambda was that we really shouldn't have anonymous namespaces in header files. Since you moved things to the .cc file, it's all good.
2 participants
In #2139 we had moved the anisotropic viscosity assemblers out of the main code and into the tests, since there was no realistic application for them at the time. #3066 then introduced a cookbook that duplicated the existing code from the test in a plugin. #6615 now needs the same code again, so I think it is time to move the code back into the main program and unify the code. The assemblers are separated from the normal Stokes assemblers, because they change things in the innermost loop of the Stokes assembler and are significantly more expensive than the isotropic viscosity ones (that was the motivation in #2139 to split them from the normal assemblers).
To make sure nothing was changed between the cookbook and the test assemblers, I first created a test for the cookbook (meanwhile fixing a bug in the cookbook prm), and then unified the assembler code. The results of the test remained unchanged.
The anisotropic viscosity assemblers are still not accessible from a simple parameter file change, they need a plugin to replace the isotropic assemblers with the anisotropic ones. But at least we now dont need to duplicate the code anymore.
I have now also unified the shear heating plugins into a new heating plugin.
@Wang-yijun this will simplify your work in #6615